Bioclimatic Zoning for Sheep Farming through Geostatistical Modeling in the State of Pernambuco, Brazil

Simple Summary Heat stress (HS) is a complex phenomenon that triggers multiple animal response mechanisms that negatively impact livestock welfare and production. Thermal comfort is therefore an important subject for limiting performance loss and other adverse effects of heat stress on animal physiology in different production systems; furthermore, it is becoming increasingly important in light of recent climate change scenarios. The purpose of this study was to point out different areas of Pernambuco state that are likely to be best suited to different sheep breeds. The study identified two dairy breeds (East Friesian and Lacaune) that have good potential to be farmed in specific areas of Pernambuco state. The thermal comfort indices presented in Pernambuco were favorable for the main meat-producing breeds. Abstract The Intergovernmental Panel on Climate Change (IPCC) has pointed out the high vulnerability of developing countries to climate change, which is expected to impact food and income security. Sheep farming is one of the main animal productions among the families located in the most vulnerable regions of the semiarid region of Pernambuco state, a Brazilian territory known for its high temperatures, low relative humidity, and high net solar radiation. Therefore, the objective of this study was to identify different regions of Pernambuco that may be more suitable for different breeds of sheep, based on non-parametric statistics and kriging maps of the temperature and humidity index (THI). THI values were determined based on mean annual temperature and wind speed extracted from the TerraClimate remote sensing database. Pernambuco state presented THI values ranging from 66 to 79, with the hair breeds having a high potential for exploitation in almost all territories, including the main meat-producing breeds. The East Friesian breed, a high milk producer, would be well suited to the Agreste mesoregion, a territory that, like the Pajeú and Moxotó microregions, also proved favorable for the introduction of three wool breeds (Suffolk, Poll Dorset, and Texel) known as major meat producers. The kriging maps of the THI values successfully allowed the identification of strategic development regions of Pernambuco state with high potential for sheep breeding.


Introduction
Heat stress is the result of a combination of climatic variables, including wind speed, high temperature, air humidity, and high solar radiation, that adversely impact animal welfare and productivity [1]. An animal is considered stressed when it needs to alter its a series of physiological and behavioral variables. Mendes et al. [21] used geostatistics to perform a bioclimatic zoning for the Dorper breed in Pernambuco state. However, there are no data or scientific studies involving the other breeds that are relevant for the economic development and food security of the population present mainly in the semiarid region of the state. Thus, to fill the knowledge gaps in this existing area, this study aimed to identify different regions of Pernambuco state that are more suitable for several sheep breeds based on non-parametric statistics and kriging maps of annual THI values. Such maps are strategic for development plans for the region in the near future, enhancing adaptation to climate change.

Characterization of the Study Area
This study was carried out in Pernambuco state-inserted in the NEB region, which covers an area of 98,312 km 2 , situated between the parallels of 7 • 15 45 and 9 • 28 18 S and the meridians of 34 • 48 35 and 41 • 19 54 W. According to the Köppen-Geiger climate classification, the region's climate divides into the following six classes: Am (humid tropical climate); Aw (tropical savanna climate with a dry winter season); BWh (hot arid climate); BSh (hot semiarid tropical climate with a defined dry season); Csa (hot summer Mediterranean climate); and Csb (cool summer Mediterranean climate) ( Figure 1) [35].  Throughout most of the year, the semiarid region of Pernambuco presents high temperatures, with an average annual minimum and maximum of 24 and 29 • C, respectively, reaching 35 to 38 • C during the hottest hours. The intense solar radiation in the region accumulates an annual energy of 2.2 MWh·m −2 with the average daily solar irradiation higher than 5.0 kWh·m −2 ·day −1 [36].

Climatological and Geospatial Data
The data used in this study came from the TerraClimate database, which comprises water balance and climate data, with a monthly frequency and spatial resolution of approximately 4 km (1/24 • ) and time series from 1958 to the present time [37]. The TerraClimate dataset is divided into primary and secondary climate variables. The primary climate variables are maximum temperature, minimum temperature, vapor pressure, total precipitation, downwelling surface shortwave radiation, and wind speed. On the other hand, the secondary climate variables are reference evapotranspiration (ASCE Penman-Monteith standardized model), runoff, actual evapotranspiration, climatic water deficiency, soil moisture, equivalent water snow, Palmer Drought Severity Index (PDSI), and vapor pressure deficit.
For the present study, maximum and minimum air temperature (Tmax and Tmin, • C) and mean wind speed (Ws, m.s −1 ) were used, totaling 4653 observations (sampling points) ( Figure 1) distributed throughout Pernambuco state. Data were obtained from the Climate Engine platform (https://climateengine.com/, accessed on 12 September 2022), the image processing platform from TerraClimate, as well as the point data georeferencing. Subsequently, the mean air temperature (Tair, • C) was obtained, and then the annual mean temperature and humidity index from 2010 to 2021 was estimated, as well as the fractional behavior of this index in the study region [38][39][40][41][42][43], according to Equation (1) established by [40]: THI = (6.3952 +0.08964 Tair + 0.01018 Ws) 2 (1)

Statistical Analysis
The spatiotemporal data were submitted to descriptive statistical analysis to obtain the mean, median, minimum, maximum, standard deviation (SD) and coefficients of variation (CV, %), and asymmetry and kurtosis. The percentage value of the CV was categorized as low (CV < 12%), medium (CV = 12-24%), or high (CV > 24%) [44]. Subsequently, the Kolmogorov-Smirnov (KS) normality test was also applied, using a significance level of p ≤ 0.01.

Geostatistical Analysis
Kriging is a regression method used to interpolate data that takes into account the spatial autocorrelation characteristics of regionalized variables, using a mean structure and a Gaussian stochastic process, assuming that points close together in space tend to have more similar values than points farther apart [45][46][47]. By taking into account the existence of spatial continuity, it allows data obtained by sampling certain points to be used for the estimation of points where the value of the variable is unknown.
To investigate the spatial structure of variation, geostatistical analysis based on classic semivariances was adopted, according to Equation (2), which estimates the degree of spatial dependence between the pairs of observations. The magnitude of the semivariance between two points depends on the distance between them, implying smaller semivariances for smaller distances and larger semivariances for larger distances [48,49]. The plot of semivariance as a function of distance to a point is called a semivariogram.
Semivariance calculation, based on Equation (2), semivariogram function model fitting, and cross-validation were performed using the geostatistical software ArcMap 10.5 from Environmental Systems Research Institute (ESRI).
where γ(h) is the experimental semivariance estimator, obtained by the sampled values Z(X i ), Z(X i + h); N(h) is the number of measured value pairs separated by the vector or lag distance; h is the distance between sample pairs; and Z(X i ) and Z(X i + h) are the values of the i-th observation of the regionalized variable, collected at points X i and X i + h (i = 1, . . . , n), separated by the h vector.
The raw data were imported to ArcMap, and experimental semivariograms were calculated. Three variogram models (i.e., experimental, Gaussian, and spherical) were fitted to the experimental semivariogram [50].
Spherical Model: Exponential Model: Gaussian Model: where γ(h) is the experimental semi-variance estimator; C 0 + C is the sill; C 0 is the nugget effect; C is the variance dispersion; h is the distance between sample pairs; and a is the range (m).
As an alternative method for evaluating the model's accuracy, deviations in the estimates from the measured data were compared by cross-validation [50,51]. This comparison of performance between the models was carried out using the following statistics: mean absolute error (MAE), mean error (ME), mean square error (MSE), average standardized error (ASE), root mean square error (RMSE), and root mean square standardized error (RMSSE). The five error statistics of predictions were applied to the cross-validation analysis. The equations are as follows [52]: whereẐ(x i ) is the predicted value, Z(x i ) is the observed value, N is the number of values, and σ i is the standard error for location x i . The degree of spatial dependence (DSD), when less than 25%, was considered strong. Between 25 and 75%, it was considered moderate, and once greater than 75% it was considered weak [53].

Boxplot Analysis
To identify outliers and some statistical properties, a boxplot analysis was performed in the 12-year series THI data used in this study, constituting the summary of 5 numbers (Figure 2), that is, the minimum and maximum values and 3 percentiles (median and interquartile range). The resultant boxplots show that there is a significant pattern found in the dispersion of the index measured in this study. The data were homogeneous [54], and no outlier was found. The median was similar for all values, without a large discrepancy between the years, as can also be seen in Table 1. The result of the boxplot analysis indicates little asymmetry in the data studied [55], also shown by the position of the median line. The amplitude was similar for all the years studied, except for 2018, a year that was characterized by a great drought in the NEB [56,57], which showed a smaller amplitude when compared to the other years studied. Table 1 presents descriptive statistics of annual THI values for the studied years, based on a general average. The CV was considered low (< 24%) for all the years evaluated [44], indicating a high homogeneity of the data. As the mean and median values were similar, a behavior verified for all years evaluated in the present study (Table 1), the data were assumed as presenting normality, in accordance with Silva et al. [43,58].  The resultant boxplots show that there is a significant pattern found in the dispersion of the index measured in this study. The data were homogeneous [54], and no outlier was found. The median was similar for all values, without a large discrepancy between the years, as can also be seen in Table 1. The result of the boxplot analysis indicates little asymmetry in the data studied [55], also shown by the position of the median line. The amplitude was similar for all the years studied, except for 2018, a year that was characterized by a great drought in the NEB [56,57], which showed a smaller amplitude when compared to the other years studied. Table 1 presents descriptive statistics of annual THI values for the studied years, based on a general average. The CV was considered low (<24%) for all the years evaluated [44], indicating a high homogeneity of the data. As the mean and median values were similar, a behavior verified for all years evaluated in the present study (Table 1), the data were assumed as presenting normality, in accordance with Silva et al. [43,58].
The criteria that the RMSSE should be close to 1 and that the RMSE and ASE values should be close to each other led to the Gaussian model being chosen and then being adopted for the kriging maps ( Figure 3). In Table 2, we can see that for all years the Gaussian model presented RMSSE greater than 0.91, while for the other models the value was between 0.52 and 0.72. The determination coefficient (R 2 ) showed an adequate fit, with values greater than 0.9, for all the years evaluated, as shown in Figure 3. The parameters (i.e., nugget effect, sill, and range) as well as the degree of spatial dependence of the geostatistical model used can be seen in Table 3.

THI Kriging Maps
Using the validated and established semivariogram models, maps of the spatiotemporal distribution of annual THI for Pernambuco were processed ( Figure 3) by kriging. To facilitate the bioclimatic zoning, the map of Pernambuco state was divided according to the physiographic zones of its own mesoregions: I-Metropolitana; II-Zona da Mata; III-Agreste; IV-Sertão; and V-São Francisco.
The maximum variation in THI across the state ranged from 66 to 79, where the higher the value, the more stressful the environment [60]. However, the threshold value, at which point the animal starts to be stressed, depends on the species and even the breed [19,61,62]. The São Francisco mesoregion presented the highest THI values in the state, while Agreste was characterized by the lowest values. The Metropolitana mesoregion and a large part of Zona da Mata are coastal regions; therefore, they are subject to high temperatures associated with high relative humidity [57,58,63], a characteristic that may have contributed to the high THI values in these regions, especially in the eastern part of the Metropolitana area. These results corroborate with the study carried out by Mendes et al. [21], which identified that for the coldest (June) and the hottest (January) months, the Zona da Mata and Metropolitana regions had one of the highest temperature values and the highest humidity values for the state. In addition, these two mesoregions, Metropolitana and Zona da Mata, have the highest rainfall rates in Pernambuco state, with annual rainfall depths higher than 1200 and 900 mm, respectively [58].
In order to identify which regions of Pernambuco state would be most suitable for different purebred sheep, the critical THI limits per breed for environmental controldetermined by McManus et al. [19] in an assessment of the distribution of sheep herds in Brazil and their relationship to climatic and environmental factors (Table 4)-were adopted.
We also calculated the frequency distribution of the THI values of the presented mesoregions by their cumulative distribution functions, presented in quadrennials, as shown in Figure 4. According to Figure 3, there is high variation in THI values within the same mesoregion, especially in the larger ones such as the Sertão, São Francisco, and Agreste mesoregions. To enhance the discussion of the present study, Figure 5 shows Pernambuco state already divided into its microregions. The state is subdivided into nineteen microregions [64]; however, one of the microregions is Fernando de Noronha Archipelago, known for its high touristic value, which was not taken into account in the present study.  According to Figure 3, there is high variation in THI values within the same mesoregion, especially in the larger ones such as the Sertão, São Francisco, and Agreste mesoregions. To enhance the discussion of the present study, Figure 5 shows Pernambuco state already divided into its microregions. The state is subdivided into nineteen microregions [64]; however, one of the microregions is Fernando de Noronha Archipelago, known for its high touristic value, which was not taken into account in the mesoregion, especially in the larger ones such as the Sertão, São Francisco, and Agreste mesoregions. To enhance the discussion of the present study, Figure 5 shows Pernambuco state already divided into its microregions. The state is subdivided into nineteen microregions [64]; however, one of the microregions is Fernando de Noronha Archipelago, known for its high touristic value, which was not taken into account in the present study.

Hair x Wool Breeds
From Table 4, it can be seen that, except for the Bergamácia breed, the hair breeds have a higher tolerance to heat stress compared to wool breeds, based on THI values alone. This may be due to the fact that the hair structure protects the skin against direct solar radiation while promoting convection and heat loss through evaporation [65]. On the other hand, the wool cover makes water evaporation from the body more difficult, thus reducing heat loss through transpiration, although it also offers protection against direct solar radiation [66]. Therefore, the thermoregulatory process tends to occur more slowly in wool sheep [67]. However, it is not only the type of covering that influences the tolerance to heat stress of the animal. The adaptation to aggressive environments, with important physiological and structural changes, such as energy metabolism and body size, also determines the degree of adaptation of the sheep species to this type of stress [25,68].
According to the THI maps (see Figure 3), the hair breeds Rabo Largo, Morada Nova, Somali, Cariri, Santa Inês, and Dorper, as well as the wool breed Bergamácia, can be used in any mesoregion of the state. We can see that for the three periods used for cumulative distribution analysis, the mesoregions Agreste, Sertão, Zona da Mata, and Metropolitana presented 100% of THI values below 77, and São Francisco about 95%. From Figure 3, we can identify this region of THI > 77 as the westernmost portion of the Petrolina microregion, a territory where Santa Inês and Dorper sheep might be working close to their established THI limits. The White Dorper species has the lowest THI limit among hair sheep, so the extreme region mentioned might not be suitable for this breed.
The Rabo Largo sheep are known for their ability to walk long distances and to cope with harsh environmental conditions such as long periods of drought and high temperatures [69]. Although herds of this breed can be found in the São Francisco mesoregion [19], there is high potential for further exploitation of this breed in Pernambuco state, since it is usually chosen for the production of good quality meat which is closely related to the nutritional habits of indigenous human populations in other arid and semiarid parts of the world [70,71].
The Somali breed is adapted to a dry climate and scarce food supply and is restricted to the NEB. However, the population of Somali sheep mainly consists of small herds belonging to research institutes and a few herds belonging to breeders [72,73]. There is not much information about the Somali sheep in Brazil; most studies with this breed are based on crossbreeding to produce more meaty animals [74]. Furthermore, Bergamácia sheep have been considered robust, with lower maintenance requirements and easy handling; however, because of their smaller size, they are often considered to be less productive [17]. In the NEB, this breed was used mainly for intercurrent crossing with indigenous breeds [75].
Due to the replacement of local genetic groups by improved exotic breeds, the Bergamácia breed is practically no longer used by breeders [72].
McManus et al. [19] showed that wool breeds are mainly limited to the South and Southeast of Brazil (SEB), regions that have a more temperate climate and less aggressive environments. From Figure 3, it can be seen that the central area of the mesoregion Agreste (including the microregions of Garanhuns, Brejo, and the southern part of Ipojuca municipalities), as well as a small part of the extreme west of it, would be suitable for all wool breeds. It can be identified from Figure 4 that in the period from 2010 to 2017, the mesoregion of Agreste presented about 85% of THI values below 72, while for the 2018-2020 period, about 75%. From Figure 3, it is verified that the microregions of Vale do Ipanema and Médio Capibaribe presented THI > 72, extending to the region of Alto Capibaribe in the last quadrennium studied. These regions might not be suitable for Border Leicester, Lacaune, Merino, Corriedale, Ilê de France, Karakul, Creole, and Hampshire breeds. The eastern part of the mesoregion of Sertão (e.g., Pajeú and Moxotó) may also be suitable for the SAMM, Suffolk, Poll Dorset, East Friesian, and Texel breeds. However, this region is notorious for its extreme drought events [57,76], where less suitable animals may not properly develop.
The SAMM breed is an efficient feed converter and does extremely well in feedlot and pasture systems because of its ability to utilize low-quality roughage [77], presenting favorable growth and meat production attributes even under intensive rearing conditions [78]. It is one of the most common feedlot breeds throughout South Africa [79], being one of the most heat-adapted wool breeds, with potential for breeding in the recommended regions of Pernambuco state.
In addition to having low tolerance to heat stress, based on THI values (Table 4), the Border Leicester, Merino, Corriedale, Karakul, and Creole wool breeds are not recognized for being major producers of meat or milk [80]. For this reason, they were not considered feasible breeds for further exploitation in the state. McManus et al. [20] showed that, based on respiration rate, the Merino, Corriedale, and Creole breeds would usually be under stress in the semiarid region of Pernambuco.

Main Meat Production Breeds
The most popular breeds used for meat production include the Morada Nova, Santa Inês, Dorper, Suffolk, Poll Dorset, Texel, Ilê de France, and Hampshire breeds [81][82][83][84][85][86][87][88]. As mentioned, the Morada Nova, Santa Inês, and Dorper varieties can be raised in any mesoregion of the state, based on the THI maps ( Figure 3). However, in the Petrolina microregion and most of the southwestern part of the Sertão mesoregion (south of the Araripina and southeast of the Salgueiro microregions), they may be at risk of being under heat stress on an average day, especially the Dorper sheep.
Morada Nova sheep represent one of the main sheep breeds in the NEB, since they are small in size and resistant to semiarid conditions, and this breed is a main source of protein for rural populations and small farmer holdings [89]. This breed uses a thermal storage mechanism to retain heat from the hottest times of the day, releasing it during the night or early morning, and they use panting as an immediate response to environmental heat stress, with sweating being a secondary response mechanism [90]. Santa Inês is the largest breed in the country, being reared in the Northeast, Midwest, and Southeast of the country [19,74], with Pernambuco state having one of the largest herds in Brazil [91]. It is the result of crossbreeding between Bergamácia, Morada Nova, and Somali, as well as other breeds, with no definitive standard, and is a dominant breed for meat production [92]. Based on the THI limit values (Table 2), the Santa Inês breed is not among the most heat resistant, in sixth place. However, Titto et al. [67] showed that after sun exposure, the increase in rectal temperature and respiration rate of Santa Inês sheep was significantly lower than the Morada Nova breed sheep, one of the most tolerant when taking into account only THI values.
Despite being an exotic breed in Brazil, Dorper sheep performed similarly to localized sheep (Morada Nova and Santa Inês) when exposed to heat stress [93]. In addition, Dorper sheep have coat characteristics that favor less thermal insulation and greater resistance to solar radiation [94]. Costa et al. [95] observed that in an environment exposed to direct radiation, Dorper sheep showed an increase in the area occupied by sebaceous glands in the dermis.
Temperate climate breeds such as Suffolk, Texel, and Ile de France show far worse physiological responses to heat stress compared to breeds adapted to the semiarid tropics such as Santa Inês and Morada Nova, representing extremely low heat tolerance [67,74,96]. From Figure 3, it can be seen that only the central region of the Agreste mesoregion (including the microregions of Garanhuns, Brejo, and the southern part of Ipojuca) and the small westernmost part of it are suitable for these breeds. However, because these breeds are less adapted to high temperatures, they may be living close to their established THI limit, as previously mentioned.
It is suggested that high-meat-production breeds are more likely to exhibit adverse effects of heat stress due to lower heat adaptation capacity as compared to resistant breeds that are more adapted to heat [88].

Main Milk Production Breeds
East Friesian and Lacaune are the two most common dairy breeds worldwide [97,98], both being wool breeds, with the East Friesian being slightly more heat resistant, based on THI values (Table 4). According to Figure 3, the central region of the Agreste mesoregion (including the microregions of Garanhuns, Brejo, and the southern part of Ipojuca), as well as a small part of the extreme west of it (northwest of the microregion of Vale do Ipanema), may be well suited to both of these breeds. The extreme north of the Pajeú microregion may also be suitable for East Friesian sheep.

Implications of the Study
The sheep industry in Brazil is growing, and research is needed to help breeders effectively increase investment in the sector. This study identifies the environmental conditions that are favorable for sheep production in Pernambuco, showing that breeds of interest that are not explored within the state can be inserted, allowing the local economy to take advantage of their good production levels, and also serving as a warning, and guiding producers so that breeding does not happen randomly, disrespecting the specific requirements of each sheep breed.
The association between the environment and production is dependent on technologies and investments, but these resources are not always available in Brazil, especially in the semiarid part of the NEB, and the absence of policy and guidance can negatively affect production, exposing the need for a management plan for sheep farming in Pernambuco that takes into consideration the environmental characteristics of the state, as it comprises a massive part of the entire Brazilian sheep herd.
Obviously, THI is not the only factor that limits the use of a breed in a given situation; however, comprehension of the spatial distribution of breeds, which is highly correlated with environmental control, can assist in the creation of environmental descriptors, classifying breeds for conservation. This is a first attempt to zone animal production by breed in Brazil. The development of a typology based on the methodology employed in this study can promote the creation of policies intended to increase the effectiveness of sheep farming, as the ability to adapt to climate change can determine whether producers, states, and countries will increase, maintain, or decrease their production levels or market share over the following years.

Conclusions
This study identifies the potential and thermal comfort conditions in different regions of Pernambuco state for different breeds of sheep found in Brazil, based on THI parameters and using geostatistical methods. We highlight that the THI values decreased from the Metropolitana mesoregion towards Agreste and increased towards Sertão and São Francisco.
Pernambuco has the proper conditions for the production of hair sheep, including three major meat-producing breeds (Morada Nova, Santa Inês, and Dorper), for almost all of its territory, except for the most southwestern part of the Petrolina microregion, a territory marked by high THI values during the entire studied period.
When it comes to wool species, the mesoregion of Agreste presents the greatest potential for breeding, especially the more adapted and better wool breeds for meat production such as Suffolk, Poll Dorset, and Texel, because almost the entire region showed THI values below 74. However, the less adapted wool breeds, such as Ilê de France and Hampshire, which are also good meat-production animals, may still face thermal stress in the microregions of Vale do Ipanema and Médio Capibaribe. The East Friesian sheep, one of the main dairy breeds, can also be well adapted in the mesoregion of Agreste, including the region known as the Dairy Basin of Pernambuco state. Data Availability Statement: Publicly available datasets were analyzed in this study. These data can be found here: https://climateengine.com.